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A fast computer program ior predicting nonequilibrium rocket plume 
properties is described. The analytical model assumes parallel turbuient (vr 
iaminar) mixing between concentric chemically reacting streams and canaiso 
be used for studying chemical lasers and re-entry wakes. The equations for 
iree shear layer mixing with nonequilibrium chemistry are solved via a 
mixed implicit/explicit finite difference scheme which efficiently predicts 
flow preperties and composition, even when many chemical reactions are 
near equilibrium. The stability problems inherent in fully explicit fimite 
difference schemes are shown to be eliminated, and stable integration step 
sizes are shown to be increased by up to 4 orders of magnitude. Computer 
run times for typical afterburning rocket plume calculations are shown to be 
decreased by more than one order of magnitude over the original (fully 
explicit) AeroChemaxisymmetric mixing with nonequilibrium chemistry pro- 
garam. Both the analysis and computer program write-up are presented, 
including a sample calculation and a FORTRAN listing. : 
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ABSTRACT 


A fast computer program for predicting nonequilibrium rocket plume 
properties is described. The analytical model assumes parallel turbulent (or 
laminar) rnixing between concentric chemically reacting streams and canalso 
be used for studying chemical lasers and re-entry wakes. The equations for 
free shear layer mixing with nonequilibrium chemistry are solved via a 
mixed implicit/explicit finite difference scheme which efficiently predicts 
flow properties and composition, even when many chemical reactions are 
near equilibrium. The stability problems inherent in fully explicit finite 
difference schemes are shown to be eliminated, and stable integration step 
sizes are shown to be increased by up to 4 orders of magnitude. Computer 
run times for typical zfterburning rocket plume calculations are shown to be 
decreased by more than one order of magnitude over the original (fully 
explicit) AeroChemaxisymmetric mixing with nonequilibrium chemistry pro- 
gram. Both the analysis and computer program write-up are presented, 
including a sample calculation and a FORTRAN listing. 
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NOMENCLATURE 


defined by Eq. (15b) 

defined by Eq. (33) 

attenuation per unit length 

constant in expression for k, (see Eq. (23)) 
defined under Eq. (25b) 

activation energy {see Eq. (23) 


specific heat of rnixture (2 X; cy ) 
; : 
= 1 


specific heat of ith species 

electronic charre 

defined as X,/W (= Y;/W;) 

factor used for the external control of integration step size 


Gibbs free energy of ith species at standard stace (1 atm) 


change in standard Gibbs free energy for a reaction, ) (g;) products 


i 
= y (g;) reactznts 
i 


enthalpy of mixture 
enthalpy of ith species 


heat of formation of ith species at T = 298K 
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forward rate coefficient 
eddy viscosity coefficient (see Section V]) 


eddy viscosity coefficient for Donaldson/Gray model (see Eq. (32)) 


equilibrium constant 

Lewis number (laminar or turbulent) 

electron mass 

Mach number at half radius, defined under Eq. (32) 


temperature exponent in reaction rate equation (Eq. (23)) 
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electron density 


arbitrary dependent variable (see Eqs. (13) - (15)) 


iv 


electrou-neutral collision cross section 


Pp static pressure 

Pe Prandtl number (laminar or turbulent) b 

r coordinate normal to jet centerline fi 

ry inner mixing zone radius f 

se jet radius B 
ri defined under Eq. (28b) : 
r2 defined under Eq. (27) af ; 
R universal gas constant : 
T static temperature ~ 
u x component of velocity 4 
v r component of veclocity 

Me electron velocity : ; 
Ww: molar rate of production of ith species ‘ : 


w molar rate of production from jth reaction 
W molecular weight of mixture 2 F;) = 
W; molecular weight of ith species 
x coordinate parallel to jet centerline 
Ax +1 7 *n 
mole fraction of ith species 


mass fraction or ith species 


constant for external control of eddy viscosity (see Eqs. (25) - {31)) 
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yj 
¥m +17 Ve 
eddy diffusivity for turbulent flow; defined as It/p 
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n defined by Eq. (25) eA : 
p viscosity (or eddy viscosity for turbulent flow) ial 

Ve electron-neutral collision frequency : 
p density — * 
o electrical conductivity ; 


v stream function 


w Signal frequency 
i 
: ' 
Subscripts t 
aaah reece st i 
e evaluated at edge of mixing layer (free stream) 
27 electrons 
i ith species } 
} 
j value at nozzle (jet) exit ; 
m W index in grid network 
n x index in grid network | 
fe) evaluated at axis of symmetry, r = 0 
Miscellaneous 
i | 
| | absolute value a 


(ae partial derivative with respect to 8; Y being held constant 


» summation over i species 
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Continuing Solution 
Electron Density and Collision Frequency 


A FAST COMPUTER PROGRAM FOR 
NONEQUILIBRIUM ROCKET PLUME 
PREDICTIONS - Fortran Listing 
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SAMPLE OUTPUT 
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IMPLICIT/EXPLICIT AND FULLY 
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I. INTRODUCTION 


The need to accurately determine electromagnetic wave/ plume 
interactions!»* has motivated the development of several programs*~* to pre- 
dict the electrical properties of turbulent rocket exhaust plumes.* Each pro- 
gram uses the same gas dynamic model {parallel mixing between two concen- 
tric streams), but only the AeroChem program of Mikatarian and Pergament 
accounts for nonequilibrium chemistry effects; tte others?~> assume local 
‘thermochemical equilibrium to prevail. Attempts to account for transverse 
radar attenuation data (taken with focused microwaves) under simulated alti- 
tude conditions using equilibrium chemistry programs generally fail. Irideed, 
it has been demonstrated by Pergament and Jensen’ that finite-rate chemical 
kinetics must be incorporated into rocket plume calculations to accurately 
predict plume temperaiures, electrical properties, etc. 


Es st 


Conceptually the numerical solution of the equations describing 
axisymmetric mixing with nonequilibrium chemistry presents few difficulties. 
A fully explicit finite difference scheme* yields an accurate solution to the 
problem. Unfortunately, however, on many occasions the stability require- 
ments associated with the explicit integration scheme are found te so severely 
limit maximum step sizes (e.g. on the order of 1075 to 107! £t) that computer 
solutions become economically prohibitive.f This "stiffness" in the equa- 
tions,” !° (which results in small step sizes) is found to be restricted to non- 

, equilibrium flows where one or more of the chemical reactions is at or near- 
equilibrium, a situation typical of relatively low altitude (< = 70 kft) afterburn- 
ing plumes. 


Ree a SARL CLLR ONRE COROA ET TUNA SIT ONAN Fl HOES MAREN 6 


In this computer program we utilize a new numerical technique to 
solve the partial differential equations (in finite-difference form) describing 
tarbulent (or laminar) shear flows with nonequilibrium chemistry: Implicit 
differences are used in the solution of the species conservation ey. ationst 


*References 3 and 5 are unclassified descriptions of the solution techniques 


used in the Naval Weapons Center and Lockheed Propulsion Company codes. 
Further information on these codes may be found i- the classified literature. 


at 
TFor frozen (non-reactive) flow the governing equations can be integrated 
quite rapidly using explicit differences. 


SS SRET (ck Ianemenenesreet, wom nena ater enema ener a a, womens earner = 
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tThis appears to be the first program to utilize implicit differences for free 
shear layers with nonequilibrium chemistry, although such schemes have 
been applied to one-dimensional nonequilibrium nozzle flows.!!)!2 
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and explicit dificrences are used for the momentum zhd enerey equ2tions. 
This mixeé implicit/ explicit scheme eliminates the instability problems 
which characterize fully explicit schemes and allews-integration step sizes 
tobe increased by orders of magnitude without Secrificine ‘ ‘@ccuracy. Those’ 
rocket exhaust plums predictions. which could not be made with the original 
AeroChem program * (which uses the iully explicit sckeme) beczuse ihe stable 
stép size apzroached 1o77° ft, are quite easily hardled by the new programs. 
Most importantly, for the more typical czse in which average step sizes - 
range from 1G~* to 10 ft, themew program teduced computer run times by 
at feast an order of magnitude: (A g¢etailed comparisen Detween rocket plume: 
calculations using the original tuily explicit difference scheme and the rew 


mixed implicit/, explicit scheme is given in Appendix. D.) Foy 8 


” 


The gas dy namic model 2s-iumes pareliel,mixing between the rocket 
exhaust products and surrounding air ‘either quiescert or ‘moving), and aliows 
for non- uniform initial conditions at the nozzle exit plane.* Lewis and Prandtl 
numbers are assumed to be constant, and pressure. is allowed tc vary paraliel 
to. the plume axis. Turbulent transport is described via an appropriate eddy 
viscosity model and Sutherland's law is used to calculate the Viscosity for 
laminar flow.. The program is quite general; it allows any chemidéal reaction 
mechanism {and associated rate coefficients) to be used as long as thermo- 
dynamic data aze available for all species. Thermodynamic data, taken 
directly from the JANNAF Tables, ate supa intabular form. ' : 


This' report gives the governing pexiial differential equations and; 4 
their finite difference'formulations, the various eddy viscosity models which 
may be input -to the program, a description of the input data, the output from 
a sample calculation and a complete FORTRAN listing. The computer output 
gives detailed axial and radial distributions of velocity, temperature, density 
and species mole fractions? These results are used ito calculate electron 
densities and electron-neutral collision frequencies wiiich, in turn, are used 
to determine radar attenuation PANS VEEEE to the plume and electrical 
conductivity. ~- : } 

- Although the program was written in response to'a need for fast 
predictions of relatively low altitude afterburning rocket plume electrical 
properties, it can equally well be. used for applications to missile Pe heat- 
ing and IR radiation problems. . 


*Thus it is possible to account for the effects of an initial boundary layer on 


the properties in the mixing region. ; 
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H. GOVERNING EQUATIONS 


eed 


A. Conservation Equations and Beendary Conditions 


e geeswna 
‘ 


The following eqrations Cescribe the free-shear layer turbalent or 


a laminar mixins of co-flowing axisymmetric streams undergoing chemical 
££. reactions. Fer turbulent flow ali properties are interpreted to be the mean 
: (time-averaged) values. The eddy viscosity, #, is then described by one wf the 
phezomenological expressions giver: in Section VI. 
BS es 
$73 Global Cuntinuity 
€. z i 
5 Q a f-.. os 5 g oo _ : 
—.. gx Cl = > Br sy a0 (1) 
: a 
sae Conservation of Species 
sae oF 
E te. 6x” @r x @r \pr ory “i i?) | 
= ‘ 
4 \ 
EC & Conservation of Momentum | 
ie | 
ane ou du dp 1 @ eu 
ie _ = ee + —_ — anes 
ke Pa ox * OY OF dx xr ér (»=32) (3) ! 
| Conservation of Energy 
aT aT j dp fauy 128 icp aT j 
= 5 — =e—— + — fj tm ee § — f+ 
Pcs [a3 Se | “ax * pS) “x Or jpr** | : 
(4) 
2 Le oT y = ar y th 
: W Pr or /,%p, or ~ 2, “it 
= i i 
fi State 
a _ PW 
R #3 = RT (5) 
> 3 
: 4 . { 
‘pa 


baie a6 
a geann 


The conservation equziions are sclved scbject to the following initial 
and boumdary camnditions: 


x= 8: « = efrf, Fp = F(x), TF = T) 
r= @- — = — = — = 4 #6) 


roe: aT fae ae 


Pressure is allowed to vary im the exial direction according te, 
P = co = Gx + Gt = 


es (7) 


where Cp, C;, Co azd cs are impst coefficients (Card 5, Cots. 21-69). 


B. Transformztion to Stream Function Ccordinztes 


If is convenient to transform the equations into 2 streamline coordinate 
system and utiiize the stream function, ¥, 2s the redial coordinate. The trans- 
formztion irom cartesian (x, rj} coordimates to streamline (x, “) coordinuies 
(“hich automatically satisfies global continuity, Eq. (1)) is defined by: 


« oY 
¥a- = pur (82) 
x Jf ? 
ar ae (3%) 
From Eqs. {€a) aad (8b} we obtain, 
a _ pur @ ‘ 
(2 x COU ONO Ay (22) 
Oye (2 _ pvr (A 
(5 | (a), y \av/, (9b) 


Introducing Egs. (9a) and (95) into Eqs. (2), (3), ane (4), gives: 


se ij % 7 
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Species 
% 
: A 
sib 2S (eye Si], 5 tes) i 
@x i ef wNPri oF oF 2” op mop f 
ag 
acd, om the exis of syummetry, r= ¥ = O i 
ae pte, 7*3 , *s 
— = 2z (=) at ae (163) 
hfomentamm 
rT 2a 
ga} © 22 8 Four Go] cs} 
Puedx< "8 ay | YW eH) x J 
and, on the axis of symmetry, r = ¥ =0 
ée 1 dd ae 
See | eet “ean at See ues . i 
ox pa @x ~ 2p ae (ib) 
Energy 
c —4 
at _ le 1 Vn at 2 fe owt or], 
‘P ox * p ax pu f/f, + 27° 8 of LPr | oq * 
i (12a) 
peer? [ eye e FT 
gw _ oes * Pr 3 i ow } 
i 
and, on the axis of symmetry, r =f =0 
or 1 ap Cp,aT 1 ° 
ax pp dx or (2 ay2 ss pu ne (12b) 


Cc. Fintte-Differecce Formulation 


The govwernina set of parabolic partial differential equstions, (Eqs. 
(16), (28) ana (22)3, are first rewritten im finite difference form and then 
soived esume 2 forward marching teckcriqus. The chemtistry terms, wy, im the 
species coxtincity eqeations are evaluates via implictt-differences; the diffe~ 
sion terms im the smectes conticaity eqestions amd the complete energy and 
mppmnemtz.m equations are evaicsted win exolictt-differemces. The following 


fintte-differer-ce fermalations™ are used: 


(=) = Ql a ED ree 
Sz Jnl, mm 4S= ' 
eQ 2,mil ~n,m-l waa 
aa) a 2AF tle} 
i 2 (22) F =n, mit CQ. omil me 
| oe oF “a, m (asy” 
(145) 
where 
Sec made n,m . n,ml (152) 
and 
pur? 
a= =~ (15b] 


The difference equations that resale from applyins Eqs. 423). (12) 2nd (15) 20 
Eas. (20), (285 amd (82) are: 
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Fhe Ji = ees : 
aa arate Fs, oth , [ee a, mz —(Fi}e, m : 


(26a) 
. {Le ? 
= oe?) m-2 le Ede, m-1 — (Fils, «i 
+ in, m 
and, om the axis of symmetry (m = 0), 
FJo:1,0 ~ da:1,04>/eednes 
44x /#Le;, ) ) 7q 
tavy Sr 7a0 [fds 2,9 | 
(165) 
= F3)o.0 
Momentum 
Ax H | 7 
u = 1]4u —u 
nzi,m mA Ty }. m3 n,m-4l n,m 
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and, om the axis of symmetry, m= 0 


ar 4 
Bose =- (] ).. = se [c, 27% al 
* 2,0 ane (2%) ? : 
{876) 
= “6 


Enerey 
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r om a ie (2) _ Ax (2) (« aes 
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« i YY Beacay-¥ ; 
* Thm (Peo m Cit, mn (182) 
and, on the 2xis of symmetry, m= 0 
Tous0 = GES (Baas +225 (Bao [a1 ~ Fa 
(18b) 
1 ss 
* Tho (pu), 9 2 (hw) 0 
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The species mole fractions at siation #71.m are determined from 
the species conservation equations by imearising the chemistry terms* (5. e- 
(*)231 2) ad inverting the resulting matrix. The limearizations involving 
species Fy and F. {for a two-body reaction) cr Fj, F. j and Fy (fer a three- 
body reaction) at station mt} (zil varisbles are known at station nm), 2re given 
by 


FF Doss = (FiF ps a Fel aes 5 Fa, | + Fin te pnet -{F j a | 


- (iF), = jn Fine: * Fie Foc 


(FiF GFidacs = -2 (iF Fido + FF ide {Filoss + (FiFidn (F jinz 1 


(20) 


afs 


(FF Dn (Fide 1- 


The terms underscored by 2 single line contribute to the 2lements df the 
coefficiert matrix, while the terms underscored by 2 double line contribute 
to the known column matrix on the right kand side of the matrix equation for 
the Enearized system. Thus the matrix cauatien takes the form (for N 
species), 


fa 212 243 Ary | SFr 2, | 
22) F2 Q2 
E ’ 
¥ H (21) 
: ONL ann} [Fin] [x 
a i 


¥The chemistry term does not make the energy equation "stiff"; thus when 
solving for temperature the chemistry term is treated explicitly. 
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A. Inteesration Step Size 


* Am exact ¢iability relationship governing the maximum allowable 
imfesraticn step size,la=), a. cannot, be obtzined due to the nociinearity of the 
= governing equatzioas. Instead the step-size requirement to salisiy stabitiry 

: can only be estimated. Following Von Neuman” a limit is placed 02 {4x),,9. 
oe such that al? dependent variables (u, F, Fj) always rem2in equal to or srexter 
than zere. hes the maximum: step size at each radiai srid point (m 76) as 


a eee en 


1 estailished from the species conservation equztion* is estimated 10 ke 
=: m (Aw P 
a -% oe! a 4 
(Sx) cae= — = EDL (222) 
Ler 3 .2 (Pr: i 
“Rn, miz “n, M- 


Lid 
¥ 


ané on the axis,{m = 6 


FDL (22by 


iz a a a I - 
The actual integration step size, Ax . is iaken to be 5 of the smellest value 


of (Ax) 2x 25 computed irem Eos. (22a) and (22bj. This frazticn ( 5} was 
selected on the basis of many triai calculations. Since Eqs. (222) and (22b) 
ere only approximate an additional factor, FDL, has been incoz porated into 
the program in order to maintain external control (Card 5, Cots. 11-20) on 
Ax in case smaller step sizes are required to maintain stability of the solution. 
in addition, the step size can never exceed the input print incremeut (Card 4, 
Cols, 21-30}. 


Should the computed species mole fraction at any radial point Secome 
negative (fypically, because the chemistry is "fast", and one or more reac-~ 
tions are near-equilibcium), the step size is repeatedly haived until either the 
species mole fraction becomes. positive or the step size becomes less than 


* Applying the same criteria to the momentuin and energy equations usually 
results in larger values of (Ax) ay- 


terete, LL cL SRN NAL I IIT GOIN SD LNs ATOR AR te PORE LO ACA NEAL 4 Noting: § patil intel s+ tere IRE At ant hiss itll tha 206 FeO Geen ands tek BoD SE Finnigan tenets on 
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some minimum step size.* In the latter case, the program terminates. 

After the species mole fractior: becomes positzve, the next value of Ax is 
a@g2in computed from HA) nex Thus the step size is never determined fromm 
the value of Ax needed te satisfy stabilicy for the previous step. This is a 
somewhat unigue appreach to specifying Ax for the solution to finite difference 
equations, and can save substantial computer iimie. It was adopied because, 
typical rocket piume calculeatigns show that, when using the mixed implicit/ 
explicit difference scheme, the chemistry only influences stable step sizes in 
a smal] region of the fiow. Once the program integrates through this region 
step sizes of (45), oy uszaily suffice. 


B. Edece of noxing Lever 


An additional radiai mesh point (at frce stream conditions) is added 
whenever the “next-to-the-last" radial point value of témperature or velocity, 
differs from the corresponding free stream value by more than a specified 
percentage of the free stream value. One percent has been selected for 
velocity and five percent for temperature. 


CG. Halving the Mesh 


The number of grid points carnot be aliowed to expand without bounds 
because of the limited storage capacity of the computer. Therefore, the num- 
ber of points is hatved cither when the mesh increases to twice its original 
size (Card 2, Cols. 1-5) cr the number of points exceeds 26. The computer 
prints all output at the station at which halving occurs. 


*The minimum step size is an input number {Card 5, Cols 1-10). 


ae ae 


” ee, Bias. vl Sl a a ac SRS Rae Rat E-alert Be ao ie aie i 
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z iv. CHEMICAL REACTION RATS EQUATIONS 
= ; Ten possible reaction types are included in the program: 
Reaction Type 
: {1) A +B =z GC+D 
os (2) AtBtMz= CHM ' 
em OY (3) A+B = C+D+H#tE 
y (4) A+B = CG 
= (5) A+t+M z C+D++#M 
ie = (6) A+B eC -$ D 
an {7) AtTB+M—~- C+M 
‘= (8) A +B -~C+Dt+HE 
E (9) A+B -C 
a {10) AtM ~-~ Cc +D+M 
a Reaction types (6)-(10) correspond to reaction types (1)-(5), but proceed in 


arbitrary third body. In this program, all species are assumed to have equal 
third body efficiencies; thus, in evaluating w‘J ) Fy = (W)-'. The net rates 
of production for all reactions are written below, in the form, 


(3) o ppd) — pyli)« 


% He 
k,p* F UF 
é ae CD 


3 2 
me) _ k, p FAF, 2 kp Bo 
WwW Kp WRT 


the forward direction only. In Reactions (2', (5), (7) and (10), Mis an 


ke F F_F_RT 


cpr 
(3) WS) = ket F, Fy - 2 
p 
. kp? F 
- (j) 2 f Cc 
{ = Soh Teena, 
(4) w aie Na Kp RT 
2 o 
. PFA Ke FF RT 
(5) Pe W Kp W 


* The symbols RP and RM are used on the computer output, 


a 12 


man | sei- x ee ee ee ee 3 2 = 2 .4- 3B we et eh ce ee ree net gt nn 


-(j) _ £ “AB 
(7) = ae 
(8) oN) - pt RAF 
f ATB 
(9) wi) = kpe* F Fp 
: kp? # 
(10) wi). 4.4 


To reduce round-off and truncation errors rp!) and au) are computed 
separately for each reaction, All contributions to the molar rate of production 
of a given species are then computed and added algebraically to form W;. 
The forward rate coefficient, Ke, is expressed in the form, 
-N 
kp = AT exp(B/RT) (23) . 
and Kp is determined from, 


In Ky = ~ AG/RT (24) 


The rate coefficients are divided into seven types: 


Rate Coefficient Type* 


(1) kp = A 

(2) kp = AT" 

(3) kp = AT™ 

(4) ke = aT-2 

(5) ke = A exp(B/RT) 
(6) Kr = AT exp(B/RT) 
(7) kp = AT ? 


*Rate coefficient data for typical rocket plume reactions may be found, e.g., 
in Ref, 15. 
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- V. THERMODYNAMIC DATA 


The thermodynamic properties (specific heat, Gibbs free energy and 
enthalpy) for each species are taken directly from the JANNAF Thermochemi- 
i ~ A298} : 
pe ( —) and (h; —hzog) in 
tabular form as a function of temperature (Card 11). Linear interpolation is 
used to define thermedynamic properties at the local temperature. 


cal Tables,!® and input to the program as c 


A. Turbulent Eddy Viscosity Models 


The following eddy viscosity models!”~?° 


program: 


are incorporated into the 


‘Model 1! (Ferri)!” 


Initial region, * 


3 
if 
| 

VI. SRANSPORT PROPERTIES , 


@0.00137 x | pug —Paue | (25a) - 


P= pe 


Developed region, 


b= pe 


where bi is the value of r where pu = (pou, + P eu,)/2 and K is the eddy viscos- 
ity coefficient, usually taken to be 0.025. 


*Defined as region upstream of axial position where (u, — ve (uj =u.) = 0.95. 
TMost of the models contain a numerical coefficient K which must be deter - 
mined empirically, The value K = 0.025, taken from Schlicting,!® has been 
incorporated directly into the program. This can be changed by the program 
input data via an appropriate value for the additional constant, a, Eqs. (25-31) 
(Card 4, Cols. 61-70). 


aK bi |Potg — Pelle | (25b) 


14 


Model 2 (Ting/Libby)'? 


ea 


where 


" 
mo) 
an 
iT 
‘R 
K 
y 


Fp Les— uel e (22) (2) 


r 
n= 2 ( (p_/p)r? dr? 
Yo 9 


and fi is the value of n where u =(u, + u,)/2 


Model 3 


Initial re gion, 


H = pe= @0.00137 x p, Jug — a, | 


Developed region, 


where rl is the value of r where u = (u, ~ u.)/2 


Model, : 


Initial region, 


7 
i) 

"G 

mn 
t 


= #0.00137% p, fu, — u,| 


Developed region, 


Hom Pe 


eKripe|uy = ue | 


So a ee 


ee ee tn” ee ee ee 


(26) ii 
F 
(27) | 
(28a) | | 
; 
(28b) 
| 
(29a) 
(29b) | 
H 
i 
bed 
; 
f 4 
A 
15 e 
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Developéd rexzior 


s * = ~— 7% “ys ad £ 2 * 
: ah = Pe=e@Ki rt H ty -~ 2 e(-°) “) ns (26). 


oy 


Medel 6 (Donaldsén/Grayf* . : : ar 
InitiaF region, : 2 38 ' 
* 4 uw = p e= aK (=3 —.¥;) p.| aS tel/? (31a) 


For Mi 21.2 .K = 0.0468 + M1 |[— 0.0460 Mi + 0.0256 up | 


; F (32) 


8 “ 7 * 8 @ 


: Mer he K = 0,0248 | 


where, M1 is the value of the Mach number where v = {v, +u: é)/2 (i. 2., the 
half radiis). The speed of sound at the half radius, aris expressed 1 by, 


ay 


. 
Ww) 
s 


Ss RT 


1 
; phi a ne ‘ : 
, ; oa Le, (R/ Wa) | : : G3). 


where Wi and T2 are evaluated at the half radius. In Eq. (Bla), r, is a od 
inner mising zoné radius and is defined as the walne of r where (us ue)/ 
(uj — we) = 0.95. : 


‘ : ‘ i 3 


” Developed region, 


a 
\ 
' 
i 
{ 
i 
1 2 % a 
b= pes wk rip I ug, —'u,|/2 . : (31b) : ; 


#In the program, the specification of Model 5 means that Eq. (30) wal be 
used in the initial region and Model 2 (Eq. 26) will be used in the deveioped 
region. This is important for re- starting a problem in the developed region 
for which Modei'5 was gelected to run from x = 0, In this case, Model 2 must 
be specified on Card 2,-Ccis. 11-15, } 


: 
. H 
5 : x 


16. * 


Net emarsdaanaeanth anata cibihalimtdt as eal alates auathiaied titan ein ieee Rat Dette 


B. Laminar Flow 


Sutherland’s Las is used te describe the viscosity as 2 function of 
temperature. 


= 5 
z= 9.8 107% T* A(T +2 Ih _/it-sec (33) 


VH. PLUME ELECTRICAL PROPERTIES 


arena EE NE ES RR 


Electron density, electron-neutral collision frequency, unit radar 
attenuation and electrical conductivity are computed 2t 23! radial pcints ior 
zach axial print-out stezicn. 

A. Electron Density 


Me = 0.733 (10%) X,. pT mint (35) 


where p is in atm and T in degrees K. 
B. Collision Frequency 
-i = 
= = 27 2% -1 : 
Vom 4.57 (10°) p T re X; Qe. sec (36 


where p is in atm, T in degrees K and Q,. in cm*. The electron-neutral 

collision cross sections” used in the calculations and given in the following 

table are those which characterize typical solid propeilant exhaust plumes. 

\f other species contribute to the value of ¥4, the program must be modified. : 


Ss Wie sends Sg” VTE cee a ™ af a oe SS bili coche eee pena St eae * 1S, OE NRE ee NE OE Rg a fe te oO Ne ee Cert Oe RG PWR T bg TPO Gp gee i eee ee 


> 


Qe; em* 
co 2.03 (0p Fi246 (IG) 
co, 4,7 (10x 

EC 5.9 ve 

HCh 1.85 v7 

N 3.25 (l0-jz_ 

45 (20°), + 2.9 (109 


2 
—— Te - 
¥, = 6.21 (10%) 7 sec 


Gs; Transverse Redar Attenvation 


dbfin (37) 


where my is in mi}, ¥ in sec™ and o is radf/sec.* The program then com- 
putes transverse (normzi to the axis) radar attenuation 2long a “line-of-site" 
{i.e. radzr bearn “spot size" much smaller than the electrical plume 
diameter) from, 


2 


A= 2 af ér db (38) 


where r is in inches. 


D. Electrical Condustivity 


e*n 


0 = 2.54(10%) mho/in (39) 


e ™e 


where e is 1.6 X 10-7 coulomb and m, is 9.1 X 1077! kg. 


i 
| 


*Attenuation calculations can be made for a maximum of six signal frequencies. 


Vi. PROGRAM OPERATION 


Tke prozram, as Kisted in Apperdax A, mst be ren om a CDC 
compier since a library rectine mamed SECOND (called im scbractime TiCHy 
3s wsed to sense clansed time from the start of exectiion. This roctine 
enables the continaine solctiun cards 10 be peached when the execetien tame 
exceeds the firme ingct om Caré 2, Cols. 36-20--if pinch optiors 2 or 3 are 
selected (Card 2, Cols. 31-35). The program may be ren on a differem 
machise by remorving the sebroaftice amd replacing it with the demmy roetzc: 
shown below: 


SUBROUTINE TICK (5553 
3IJJ=0 

RETURN 

END 


When this is dene the maximum execution time input on Card 2, Cols. 36-20 
will not trigzer the punch optioz. 

Input and oviput are on the standard tave units (i.e. = for READ. 6 
for WRITE and 7 for PUNCH}. 


B. Sense Switch Cortrol 


1. Sense Switch i - When Sinse Switch 1 is on, the current values of x, 


Ig, and Neg are printed on-line. 


2 Sense Switch 3 - When Sense Switch 3 is on, the case is texminated 


ae 


and ail output is printed at the current axial ztation. 


C. Pregram Optimization 


it is recommended that the program be compiled using the highest 


level of optimization. On the CDC 6600 this corresponds to the FAST OBJECT 


CODE MODE. 
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IX. PROGRAM INPUT DATA 


The inoet data cards are explained Belew. A listing of the commcter 
Program, 2 seminle ingct data sheet and 2 sample oxtpel are giwen im 
Appendices A, Band C, respectively. 


Card 
No. Colem:,:s5 Description Format 
i 1-72 Rum identification 1246 
2 £-5 Enifial mamber of grid peints* in W & 
&-26 Number of species (22 my:ximum) 5 
BE-35 ¥Yiscosity option key 5B 
3-1: Leminzar (Ee. 34) 
Q@: gp = constant 
1z Model 3 (2a. 25) 
2: Model 2 (Ea. 26) 
3: Model 3 (Eq. 28) 
: 4: Riedel 4 (Eq- 29) 
S. 5: Model 5 (Eqs. 25 and 39) 
6: Model 6 (Ea. 31) 
i 16-20 Number of reactions (46 maximum) 5 
21-25 w, output option 5 


if 0: w. for each species is not output 
lz w. for eack species is output at 
ali redial pcints whose tempera- 
ture is greater than tke kinetics 
cut-olf temperature (Card 6, Coils. 
61-70) 


| 

| 

| oe 
coerdinate (26 maximum} 
| 

| 

| 


*Recommended number, 13 


20 


Description 


No. Coleurms Format 
z 26-30 RP, RM exipct option 5 
HY O: RP, EM for each reaction is moe 

Omipes 

iz RE. BM for exch reaction 25 occt- 
ont at ail radial points whose 
temperciure is grecter than the 
Kimetics cut-oif temmeratare 
(Card 6, Cals. 61-70) 


31-35 Card pumch option 
E €; No cards pomcked 
Iz Cards with radial distributions of 
electron density anc collision 
frequency at all printout siations 
are punched 
2: Cards for continuing 2 given solu> 
ion are puncked 
3: Both the above sets of cards are 
puncked 


i 


36-40 faximum computer execution time 5 
{ rain) (EH card punch option 1, 2 or 3 
is selected, cards will be punched after 


= « al 
this amount of computer time) 


Pressuze option 5 
if O: Pressure is 2 constant 

1: Pressure is 2 functioz of axial 

distance 

Signal frequencies (MHz) - 6 maximum - E10.3 
for which transverse attenuation calcu- 
lations wiil be made (see Section Vil). 
Leave blank if no attenuation calculaiions 


desired 
51-60 


*In practice if this option is seiected, the calculations will terminate at an 
execution time 30 seconds less than this value. 
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Section Ill is to be used. 
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ee se Elem. 


: Bo 24 - . bs fat _&% 7 sere! 
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Ne. Calamms Descriptions pe Oem i 
4 B-HQ 09 -- Esta? waktoe oF x Qn) E103 

; BN-2000 0 - Firzi velne of x (4) EOS 

22-30 © Print imeremens (f0} E103 

31-4 Lewis momber Eas 

i 41-50 Brandi number E1023 

| 51-62 Nozzle (fet) redizs (&) E10.3 

“ 

i 61-78 a, factor eased to vary eddy viscesity* Elas 

(see Section ¥2) 
5 1-10 Ax ., Minimum integration step size E18.3 
(ft) T 

11-20 EDL, ciffesion step size factor? £10.3 
21-30 = cg (zt) E10.3 
31-28 cy (at /#) El9.3 
6 Pressure Coefficients m 
41-50 cz (atm/fe*) (see Section Il) £10.3 

51-62 cz fetm/it’) Ei0.3 

6 1-30 ressure at initial value of xf{atm) E10.3 

i 11-26 Temreratre at jet centerline (K) E10.3 

21-30 Temperature at edge of jct {K) (free £10.3 

i stream value) 

4Set a= 1.0 if no changes to eddy viscosity coefficient are desired. 

Zz 

fRecommended value, 1 X 107!° ft. 
{Set FDL = 1.0 if maximum integration step size criterion discussed in 


dat 


RAE? 
BA a meres 


_No- Cefearmrs Deexcriptica 


é 31-29 Velocity an jet centerline (it/sec) 


23-50 Welecity at edge of jet (iz/sec) (free 


stre2m velae) 


3 
[ais « gitom “Bye 
53-2 AW, 132.2 sec 


ii Aw= 0 (or blank) the prsaram com- 
putes AW assuming veloctiy, teémperetere 
and species mole fractions 2ré constant 
2t imétial station 


He AWis specified, cards 7.1 to 9.13 are 
read next. Hf AWis mot specified, card 
10.1 is the next card reeda with cards 7.1 
through 9.13 Seine emitted 


61-70 Kinetics cet-cfi temperature (Kixetics are 


frozen at temperateres below this value*) 


NOTE: For the purpose of illustration it is assumed that there are J3 initial 


grid points in the Y cooréinate. 


71 1-10 741), Temperature (K) at jet centerline, 
=+=0 


7.2 1-le = T(8) 
51-60 7(i3), Temperature (K) in free stream 


8.1 1-10 = xf{1), Velocity (ft/sec) at jet centerlinc, 
r=0 


1-70 ~—_u(7) 
8.2 1-10 = u(8) 


51-66 —-u(13) 


*If blank or zero, 400K is used. 


SEER AS teats TTT ETD 
cathe! uF! Liha ae ere = 


Format 


E10.3 
£10.3 


E10.3 


E16 3 


Ei9.5 


E10.3 


E10.3 


E10.3 


RET SPERM SD ONE An in, 


. a Snaginaeyon® ap Ang gp titre an snares Aner 


Card ; 
Ne. Colemn 


NOTE: For the perecse of illustration it is assumed that there are 3 species 


G.8.% 1-10 
ree 

94.2 1-10 
9.2.0 ae 
6170 

9.222 1-10 
9.73.1 ae 
6370 

9.13.2 1-106 


NOTE: Card Type 10 


10.1.1 1-10 

61-70 

10.1.2 1-10 
24 


5 
5 
‘i ‘ 

Shae BS Statement 
, rete Sal a, Ae 


Description _Format 


Mole fraction of Ist species at jer center- E10.5 
line, + =9 


Mole fraction of 7th species 2t jet center- E19.3 
Hine, r=0 


Mole faction of Sth species at jet center- E10.3 
lize, r =O 

Mole fraction of Ist species at next grid E19.3 
point 


Mole fraction of 7th species 2t next grid 
Doint 


Mole fraction of &th species cf next grid E16.3 
point 
Mote fraction of lst species at 13 grid F10.3 


point (free stream) 


Mole fraction of 7th species at 13th grid 
point (free stream) 


Mole fraction of 8th species at 13th grid E10.3 
point (free stream) 


is not required if Card Types 7, 8 and 9 are input. 


Mole fraction of Ist species at jet center- E10.3 
line, r =0 


Moie fraction of 7th species at jet center- 
line, r = 0 


Mole fraction of 8th species at jet center- 
line, r =0 


Ss seca i 


aw. 


ee 


19.2.2 


NOTE: 


Card 
No. 


11.1.1 


11.1.2 


10. 


en TR A etic lie et a 


Description 


— 


Mole fraction of Ist species in free stream 


Format 


E10. 3. 


Role fraction of 7th species in free stream 


Mole fraction of 8th species in free stre2m 


The following cards contain the thermodynamic data.* The first 
card contains the species name, molecular wz2ight and heat of 


formation. 


The second and remaining cards contain the temperatyre 


and corresponding specific heat, free energy and enthalpy for that 


species. Two teniperatures and corresponding thermodynamic data 
The input table can contain up to a maxi- 
fhe fata are input exactly as 


are placed on each card. 
mum of 30 temperature points. 


presented in the JANNAF tabies.® 


Column 


1- 6 
7-16 
17-26 


1-10 
11-20 


21-30 


31-40 
41-50 
51-60 


61-70 


Description 
eee ts 


Name of firs: species 
Molecular weight 


Heat of formation, hog,; (kcal/rnole) 
First temperature point (K) 
c, (cal/mole-K) 

Pi 


g. - hage. 
i i 
- (4. ) (cal/mole K) 


h; - hag; (kcal/mole) 
Second temperature point (K) 


c,,. (cal/mole-K) 
Pi 


1 


&. ~ hag 
- (4.— ) (cal/mole-K) 


aoe Format 


Ab 
E10.3 
Ei0.3 


E10.3 
Bi0.3 


E10. 3 


E10.3 
E£10.3 
E10.3 


E10.3 


*The order of the species must be identical to the order on Card Types 9 or 


STEN can nch eas 


bette. Selagee See 


er OSE) See 
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Sy aia Be 


Column Description 


71-80 hy; - biygg. (kcal/mole) 


: Third temperature point 


. Name of second species 


The chemical reaction mechanism for a particular probem is input on 
the last set of cards, one card for each reaction. (See Section IV.) 


No particular order is required. 


a aaa Ea Re a ile tie 


Format 


E10.3 


1-6 Species A A6 
7 + sign 
8-13 Species B (or M) A6 
14 + sign 
15-20 Blank orM 6x 
21 = sign 
22-27 Species C A6 
28 + sign (if needed) 
29-34 Species D (or M) A6 
35 + sign (if needed) 
36-41 Species E (or M) A6 
42-48 Blank 
49-50 Reaction type, 1 to 10 (see Section IV) 2 
51 Rate coefficient type, 1 to 7 (see Section IV) Il 
52-59 A, Pre-exponential factor cm-molecule- E8.2 
sec units 
60-63 N, Temperature exponent F4,.1 
64-72 B, Activation energy cal/mole F9.1 


Next Reaction 


X. PROGRAM OUTPUT 


A. Description 


Sample output sheets for the input data given in Appendix B are 
shown in Appendix C. The first page contains the key program input para- 
meters, initial distributions of velocity, temperature and species mole 
fractions and the chemical reaction mechanism and rate coefficients. The 
succeeding pages contain printouts at axial stations corresponding to the print 
increment (Card 4, Cols. 21-30). Following the listing of electrical proper - 
ties (last page of output), the integration step size and corresponding axial 
position are given for each integration step up to the next print station. 


B. Units 


A mixed system of units appears on the program output for ease in 
making additional calculations and for program check-out. 


COLLISION FREQUENCY electron-neutral collision frequency, 
“1 
sec "2 


DELTA X integration step size, ft 
DENSITY g/cm? 

ELECTRICAL CONDUCTIVITY mho/in 

ELECTRON DENSITY mi"! 


ENTHALPY mixture static enthalpy, cal/g 


HALF RADIUS/R nondimensional radial distance to 
point where u = (uy + u,)/2, ft 
(viscosity option 6 only) 


INNER MIXING ZONE RADIUS/R nondimensional radial distance to 
point where (u-ug) / (a; - u,) = 0.95, 
(viscosity option 6 only) - 


MACH NO. Mach number 


MACH NUMBER AT HALF RADIUS Mach number at point where u = 
(ug + ug)/2 (viscosity option 6 only) 


MINIMUM STEP SIZE (AX) win, ft 


OSTA Bs wa TE 


¢ 
fo. g 


~ 


MIXING RATE COEFFICIENT defined by Eq. (3la).or Eq. (31b) 


i 
(viscosity option 6 only) 
MOLE FRACTIONS species mole fractions a 
NET RATE OF PRODUCTION 
(W-DOT/RHO+U) wh u, mole/g-ft E ; 
PRESS : pressure, atm Wo: 
1 
‘ lbm -1 - ’ 
PSI stream function, E 5) sec | 
PT radial grid point number f 
R nozzle radius 
REACTION J refers to reactions listed on first 
page of output 
RM negative molar rate of production for 
jth reaction (see Section IV), mole/ 
ml-sec 
RP positive molar rate of production for ( 
jth reaction (see Section IV), mole/ ; 
ml-sec 
“ SIGN. FREQ. signal frequency for attenuation 
calculations, mHz : 
TEMPERATURE K fe 
TRANSVERSE ATTENUATION db 
UNIT ATTENUATION db/in 
VELOCITY axial velocity, ft/sec 
VISCOSITY », lbm/ft-sec 
x axial distance, ft . 
X/R : nondimensional axial distance 
Y/R nondimensional radial distance from 
axis 
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C. Card Output 


1. Gontinuing Solution ~ Due to the large amount of input data required 
to continue a solution (after being terminated at a given value of x) the pro- 
gram will punch all input except the thermodynamic data and chemical reaction 
mechanism if this option is selected (Card 2, Cols. 31-35). However, a new 
final value ofx(Card4, Cols. 11-20) must be input manually if the problem was 
initially terminated by reaching the oid maximum distance. Cards will also 
be produced if maximum time (Card 2, Cols. 36+40) is exceeded, by using the 
proper option. 

2, Electron Density and Collision Frequency ~ The complete radial 
distribution of electron density and collision frequency at all printout stations 
will be punched if this option is selected (Card 2, Cols. 31-35). 
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APPENDIX D 


COMPARISON BETWEEN MIXED IMPLICIT/EXPLICIT 
AND FULLY EXPLICIT FINITE DIFFERENCE SCHEMES 
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APPENDIX D 


Comparisons between rocket exhaust plume calculations obtained 
using the present program (mixed implicit/explict scheme) and the origina! 
AeroChem program® {explicit scheme) have been made in order to (i) demon- 
strate that both sc!.emes yield identical results and (ii) determine the differ- 
ence in computer run times between the two schemes. The case investigated 
is a typicai solid propellant exhaust containing potassium as the dominant i 
aixali metal impurity. (Typical input data for such a propellant are given in . 
Appendix B.) The programs were rur on a CDC 6600 computer. Comparisons 
between centerline and radial distributions, velocity, temperature and species 
mole fractions computed via the mixed implicit/ explicit and explicit schemes 
are shown in Figs. D1-D4. The notation k,/100 means that the rate coeffi- : 
cients for K + HCi ~ KCl + Htwhich is very near equilibrium, even at a pres- 
sure of 0.1 atm) was «rbitr. rily reduced by a factor of 100. This was 
necessary i.. order to rus: the original] (explicit) AeroChem program® without 
excessive computer time. Figures D1-D4 show that the two schemes do, in 
fact. give identical results and we are therefore contident that the new pro- 
gram properly calculates exhaust plume properties. 
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The principal advantage of the mixed implicit/explicit scheme over 
the explicit scheme is that stable integration step sizes are much larger. This 
is clearly demonstrated via the (typical) step size comparison shown in Fig. D5. 
Thus even though the computational time for each integration step in the pre- 
sent program is about 4 times greater than in the original AeroChem proyram 
(primarily because of the matrix inversion subroutine in the present program) 
the much larger step sizes still give considerable reductions in computer run 
tame. 


Table D-I shows a comparison between computer run times for the 
present and original programs for a typical solid propellant plume at sea level 
and at 50 kft. For this case, the original «AAeroChem program cannot be usca 
(in a practical scase) to compute plume properties--primarily because tk 2 
reaction K + HC) - KCl] + liis very nearly in equilibrium--resulting in stable 
integration step sizes of less than about 107° ft. Thus we record an (essen- 
tially) infinite computer run time in Table D-I. However this case can be run 
using a modified version of the original AceroChem program® which keeps the 
reaction K + HCl ~ KCl + H identically in equilibriurn. Run times are noted 
to be 30 min and !7 min for the sea level and 50 kft cases, respectively. 
Usirg the present prograrn reduces these run times to 2.5 and 2.2 min--a 
very considerable reduction. If, in addition the program is compiled with the 
6600 optimized compiler, the run tirne then is reduced to 1 min for the 50 kit 


case. 
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Additional comparisons (made at AFRPL while demonstrating the 
present program), between the original and present AeroChem programs 
for other low altitude plumes ai3;0 snowed reductions in computer run times 
by factors of from 10 to 30. 
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TABLE D-I 


CDC 6600 COMPUTER RUN TIME, MINUTES 
10 ft plume; ™ = 0.25 ft 


Sea Level 50 kft 


Original AeroChem (Explicit) 


Program? « ow 
Original AeroChem (Explicit-modified to put 
Program® K +HCi- KC1+H 
in equilibrium) 30 17 
Present Program (Mixed Implicit/Explicit) 2.5 2.2 
Present Program (Mixed Implicit/Explicit using 
CDC 6600 Optimized 
Compiler) 1.0 
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FIG. DI COMPARISONS BETWEEN CENTERLINE DISTRIBUTIONS 
OF VELCCITY AND TEMPERATURE COMPUTED 
FROM MIXED IMPLICIT/EXPLIiCIT AND EXPLICIT 
DIFFERENCE TECHNIQUES 
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FIG. D2 CCMPARISONS BETWEEN RADIAL DISTRI JTIONS 
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FIG. D3 COMPARISONS BETWEEN CENTERLINE DISTRIBUTIONS 
OF Xt AND Xu COMPUTED FROM MIXED IMPLICI?/EXPLICIT 
AND EXPLICIT DIFFERENCE TECHNIQUES 
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